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Large eddy approximation of turbulent flow 
in DSC schemes 



Abstract. Large eddy simulation of turbulent flow is given a natural setting 
within the DSC framework of computational fluid dynamics. Periodic cellular 
coarse-graining prevents the nodal flow from piling up and preserves its large 
patterns. The coarsening operations are consistent with the near-field interaction 
principle of DSC and - therefore - uncomplicated at boundaries. Numerical exam- 
ples validate the approach. MSC-classes: 65C20, 65M06, 76D05 
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L'acte de la nature est simple, et la complexity qu'elle parait avoir 
construite ... n'est que l'entrecroisement sans fin des antagonismes 
qui se sont neutralises 

Henri Bergson 

1. Introduction 

Discrete schemes have to master the cursed situation that they first create 
by breaking space into pieces. Mesh cell systems, however fine, can never be 
perfect substitute for continuous space. Already in linear algorithms they 
produce artefacts, such as unphysical spurious solutions, and more harmful 
things happen in the non-linear case. 

In fluid dynamics, when the non-linear spectral transfer properties of 
the Navier-Stokes equations (the energy cascade; cf. Pope [Poj l come into 
play beyond the transition to turbulence, eddies are locally excited, down 
to very short scale; cf. KOLMOGOROV [Koj . Eddies smaller than the cell size 
cannot be properly resolved by a mesh of realistic coarseness and thus tend 
to induce local fluctuations that artificially pile up. Large eddy simulation 
(LES) aims to discard such divergences through regularizing the flow by 
means of a suitable averaging filter, and to screen its essential large patterns 
in this way. 

LES stands for the most promising line in modern turbulence modelling. 
The name has been coined by Deardorff in 1970, who first applied LES 
methods to turbulent channel flow |Dej . Important elements had yet been 
kept ready, formerly. Already Reynolds |Re) used temporal flow averages 
that as Reynolds averages play still a role in conventional turbulence models. 
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In 1922, Richardson [RiJ proposed mesh cell averages for smoothing down 
local fluctuations. Such so called box filters stand today for one line of LES. 
Spatial filters of varied type indeed classify modern LES methods. Excellent 
insight into the state of the art is gained from the recent book of Berselli 
et al. [BILj . which not only exposes the elaborate mathematical framework 
but also addresses many items of conventional turbulence modelling (the 
prominent k— e model of Launder and Spalding [LS], e.g.). 

The present paper differs from most related work in the field in that it 
sharply separates the natural and technical aspects of turbulent instability 
and strictly focusses on the latter within the computational context in hand. 
The natural part, viz. essentially (modulo discretization) the Navier-Stokes 
equations for viscous incompressible flow, are thereby taken as a physical 
given. This is in obvious contrast to widespread use: In many conventional 
turbulence models, as in some lines of LES, turbulent pile-up is controlled 
by modifying the underlying dynamic equations or even the physical con- 
stants. Such models use, for instance, Reynolds number dependent global 
parameters - usually an eddy viscosity, e.g. - or time and space filtered 
Navier-Stokes equations or any (e.g. stochastic) variant. All this is allowed, 
of course, as long as it provides numerical data in harmony with observation. 
Also, a profound change of the inner dynamics can very efficiently delimit 
turbulence and thus lead to a convergent scheme. It is yet, obviously, tuning 
the underlying physical phenomenon rather than eliminating the technical 
cause of divergences that are manifestly an artefact of discretization. (Tur- 
bulence is still enigmatic in many respects, but divergence is certainly not 
a feature of turbulent flow - observed in the continuum of nature, e.g.) 

Therefore, the central point of attack in discarding turbulent pile-up will 
be in this paper the weaknesses of discretization, while the Navier-Stokes 
equations are directly discretized as taken from physics - which not only 
combats the desaster at the origin, but also yields validated results. 

2. DSC schemes 

A mesh cell system is not just poor surrogate continuous space, but it always 
introduces its proper structure. Cellular meshes, with that we are dealing 
here, artificially impose local cell-boundary duality upon space |Helj - and 
it is this simple fact that dual scattering channel (DSC) schemes match in 
a canonical fashion. 

DSC schemes are finite volume methods of a very specific type. In tech- 
nical detail, they are extensively dealt with in [Helj . so it is sufficient here 
to recapitulate their essential features in a crash course like manner. The 
algorithm is characterized by a two-step cycle of iteration which alternately 
updates the computed fields within cells and on their interfaces. If the up- 
dating instructions are explicit, then a near- field interaction principle gives 
rise to a scattering process interpretation known as Johns ' cycle [Helj . A 
pair of vectors that represent the same field within a cell and on its surface 
essentially constitutes a scattering channel. Equivalently, scattering chan- 
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nels are sometimes defined as pairs of distributions that 'measure' the field 
within a cell and on one of its faces. The well known primal DSC scheme 
is Johns' TLM algorithm, wherein scattering channels are visualised as 
transmission lines | JoB| . In many applications, such as computational fluid 
dynamics, the TLM picture of wave propagation fails, but one can still go 
back to the far more general DSC setup. 

Within any physical interpretation , cf . |Helj , each scattering channel 
represents a pair of scalar or vector valued distributions z p — {p , Z ) and 
z n — (p ~ , Z) which evaluate a physical field Z on the surface and in the 
interior, respectively, of a mesh cell. A distribution p pertinent to a cell 
face is called a port and p ~ its nodal image, and the two are related by 
pull-back, viz. 

(1) (p~ Z) = (po a, Z) = (p,Zo a- 1 ), 

for every Z (of class C°° , e.g. ) , where a denotes the spatial translation 
cr : M 3 — ► K 3 that shifts the geometric node (i.e. the cell centre) in the 
(centre of) the respective face. 

It follows that there exists a scattering channel representation of DSC 
states: If M denotes the mesh cell system, d( the boundary of cell £ E M 
( 9£ is naturally identified with the set of ports with nodal image in £ ) , 
and L p the span of the z p (i.e. the linear range space), then each state in 
the mesh permits a unique representation in the space 

(2) p ■■= TT TT L p c xLf , 

with canonical projections 7r^'™ : P —> P£' into the port and node com- 
ponents of cell £ . Furthermore, there is a natural involutary isomorphism 
nb : P -> P, 

(3) nb : (z p , z p ~ ) ^ (z p ~, z p ) , 

called the node-boundary map, which hence maps P p onto P n and vice 
versa. (The cell index £ is omitted here and in the following without danger 
of confusion.) 

It stands to reason that the nodal images of different ports (e.g. on dif- 
ferent faces of the same cell) can represent the same physical field in a node 
and may even coincide there as distributions - just as two ports pertinent 
to neighbouring cells obviously represent the same field on a common face, 
if they connect two channels on it. In fact, the scattering channel represen- 
tation of DSC states in @ is in general highly redundant - which may be 
utilized for process parallelization [Hel] , 

Within the algorithm, the port and node components are updated at 
even and odd integer multiples, respectively, of half a timestep r and are 
usually constantly continued as step functions over the subsequent time 
intervals of length r . 

Note that existence, not necessarily explicit construction or application, 
of a scattering channel representation characterizes DSC schemes, and that 
less redundant representations are actually used in many implementations. 
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In fact, the scattering channel representation is basically a theoretical means 
for describing (and deriving) certain DSC properties - such as the following. 

A fundamental principle - closely related to the Courant-Levi stability 
criterion - is near-field interaction. It requires that every updated state of 
a node or face depends only states (along with their history) in scattering 
channels connected to the respective node or face - the latter here being for 
once identified with its adjacent face, if such exists in any neighbouring cell. 

As a consequence of near-field interaction, every DSC process allows for 
an interpretation as a multiple scattering process in the following sense. 

Let for any process z = (z p , z n )(t) incident and outgoing fields zf n 
and z" ut be recursively defined as processes in P p and P n , respectively, 
by setting zf n (t) := z^ ut (t-L) := , for t < 0, 
and for < t = mr ; m G N : 

zf n (t) := z v(t)-nboz n out {t-~), 

(4) 1 

== z n (t+ -) -nbozf n (t). 

Then, at every instant holds z p ( t ) — nb o z™ ut ( t — ~ ) + zf n (t) and 
z n (t + |) = nbo z p n (t) + z" ut (t + |) . Also, near- field interaction 
implies that every state is only a function of states incident (up to present 
time t ) on scattering channels connected to the respective node or face. 
More precisely, by induction holds: 

Theorem. There exists a pair of functions 3? and C , defined on back in 
time running sequences of incident and outgoing fields, respectively, such 
that for every cell ( € M the process = it " o z complies with 

(5) z n Q {t + I) = K((z? n (t - MT)) pe a f i /^n) 
and the port process = 7r^ o z satisfies 

z P c (t + t) = e{(zl lt (t + I- Mr ))„|a C;MeN ) . 

(6) I 

( ' | ' short-hand for 'in any of the (i or S ) cells adjacent to' ) 

Remarks 

(i) The statements immediately imply that z ™ out and z^ 1 in are themselves 
functions of states incident on connected scattering channels, since 

n~ 

ZCoutit + 2 ) = ^(( z m ( i -M T ))peac;MSN) - nbo z c P m (t) and 

*&».(*) = C((«^ t (t - ^-MT))n| p;/ , eN ) + 

- n& ° *c!out (* - 2 ^ 

(ii) 31 and C are called, respectively, the reflection and connection maps of 
the DSC algorithm, and afield excitation may be implicit in C , cf. |Hel] . 



(iii) Near field interaction implies computational stability, if the reflection 
and connection maps are a-passive , i.e. contractive in this sense |He2j . 
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3. Coarse-graining 

Mesh cell systems import artificial defects into space - which yet sometimes 
carry their remedy in themselves, at least in part: For instance, the cellular 
mesh, by separating the virtually unresolved cell interior from the coarse- 
grained cell-boundary skeleton grid, necessarily fixes - and hence delimits - 
the scale of local resolution. A clever design always takes into account the 
coarsening effect of the cellular mesh - which, on the other hand, provides 
quasi free of charge a basis for large eddy approximation in the DSC setup 
of computational fluid dynamics. In fact, already intuitively DSC schemes 
go along with kind of natural large pattern approximation that - in a sense 
now made precise - associates cell averages to some nodal quantities. 

The following definition is independent of any particular application and 
therefore given without reference to fluid dynamics. 

Definition (coarse-graining) 

For any port p G 9£ with nodal image p~ let B p denote the set of all ports 
on d( , the nodal images of which coincide with p~ as distributions. Also, 
let w r G [0,1]; r G B p , be a set of weights such that reBp w r = 1 ■ 

Then, every DSC field which in the scattering channel representation on 
each component q G B p C d( equals ( z q , z q ) with 

(7) z q ~ = V w r z r 

is named an (in Q and with weights w r ) coarse-grained field (pertinent 
to p - or also to its nodal image p~ ). 

The nodal state of every such field, coarse-grained in any cell, is hence 
a convex superposition - i.a. a weighted mean - of states that represent the 
field on the cell surface. Also, substituting a DSC field with any coarse- 
grained field that (obviously) leaves the cell boundary states unchanged is 
then called coarsening this field (in any component). 

Should such an innocent looking procedure untie the Gordian knot of 
turbulence ? Remember that our modest aim is repairing some defects of 
discretization - rather than resolving the fundamental questions of turbu- 
lence in the way mathematical physics should do that. What matters here 
is that in the DSC setup of CFD well-timed periodic coarse-graining of the 
flow efficiently prevents turbulent pile-up, all with retaining the large ed- 
dies. Note that the coarsening period should be taken as small as necessary 
to ensure stability and precision, but large enough against the time step 
(by one order of magnitude at least - as a rule of thumb) : Clearly, coarsen- 
ing inevitably also interferes with the flow dynamics. However - since the 
skeleton field on the cell boundaries is left unchanged - the perturbation 
is minimal, and in fact negligible, for suitable periods. Through removing 
only the turbulent fluctuations of very short range (below a scale given by 
cell size) which tend to artificially pile-up, coarsening regularises the flow 
to just such a mild degree that its essential large patterns are preserved. 
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4. Viscous Boussinesq-incompressible flow 

The port and node distributions of a DSC algorithm can be finite inte- 
grals - as is the case for the TLM method, where finite path integrals over 
electric and magnetic fields are evaluated in a discrete approximation to 
Maxwell's integral equations [He4-6]. In the case at hand, they are simply 
Dirac measures that pointwise evaluate the fields within the cells and on 
their surfaces, and composites of Dirac measures which approximate the 
gradients of these fields; cf. section 

In the Oberbeck-Boussinesq approximation Obb.BssJ all material fluid 
properties are assumed as constant - except fluid density which only in the 
gravitational term varies linearly with temperature. The energy equation of 
a Boussinesq-incompressible fluid of velocity u with thermal diffusivity a , 
heat source(s) q , and negligible viscous dissipation and heat capacity is 
then the convection- diffusion equation for the temperature T , cf. [GDN , 

dT 

(8) — + div(Tu) = a AT + q . 

The Navier-Stokes momentum equations for ( the k-th component of ) a 
Newtonian fluid of dynamic viscosity rj , under pressure p , and in a gravi- 
tational field of acceleration g take the form [e.g. Ptk] 

d 

(9) -^{gu k ) + div(gu k u) + gradp k = r/Au fc + gg k r. 

With g oc — const and g(T) = g 00 /3(T(t,x) — ) , wherein f3 = 
Q d Q / dT, this becomes in the OB-approximation 

(10) ^± + div (u k u) + ^^ = -lAu t + /?(T(t,x)-T 00 )g fc . 

at goo goo 

Integrating these equations over cell £ with boundary dC, and applying 
Gauss' theorem to integrals over A = div grad and div ( / u ) yields with 
time increment t the following updating instructions for nodal T and u k , 
these quantities and q averaged over the cell volume Vq 



T(t+-) = T 



(11) 
and 
(12) 



+ rr (agradT - Tu) ■ dS 



u k (t+^-) := u fe + t(J3(T - T oc )g k - gradpk ) + 

<s goo 



T f T] 

TF ( 9 rad u k - u k u ) • dS 

Vc Jdc Qoo 



In equations (jTTJ) and (fT2"l) the former updates (at time t — r/2) of the nodal 
quantities enter the right-hand sides in the first line, and the last former 
updates (at time t) of the cell face quantities enter the second line. 
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The nodal values of T and uj, are updated at the reflection step - 
after coarsening the velocity field, e.g. with weights proportional to the 
pertinent cell face areas, whenever the coarsening period is attained - while 
the port quantities that enter the surface integrals at the right-hand sides 
are updated on the connection step of the iteration cycle. The next section 
outlines how wc can proceed with that in a non-orthogonal hexahedral mesh. 



5. The non-orthogonal hexahedral cell 

The physical interpretation of a DSC algorithm associates a smoothly vary- 
ing (e.g. in time and space C°°- ) scalar or vector field Z to port and node 
states z p and z " of a mesh cell system. 

Let any hexahedral cell be given by its eight vertices. Define then edge 
vectors (u£)u=o,...,ii, node vectors (^6)^=0,1,2, and face vectors ( l /) t =o,...,5) 
using the labelling scheme of figure Q]a 

1 



(13) 



and J 



4 

(-1)' 



Z_^„ =0 C4M+.0" 

( (8 + 2t) e +(9 + 2(,, + (-l)')) e ) A 

A ((4+20 6 + (5 + 20 e ) L 



M = 0,1,2 



0. 



with all indices understood cyclic modulo 12 and A denoting the cross 
product in R 3 . At every cell face l s {0, 5} and for any given r 6 R+ the 




Fig. 1. Non-orthogonal hexahedral mesh cell. 
(a) Edge vectors, (b) Node vectors, (c) Face vectors. 



following time shifted finite differences of Z in directions ( /1 = 0, 1, 2 ) 
form a vector valued function 



(14) 



<\L-r 



if fj, = [t/2] 
if M ^ [^2] 



( [ x ] denotes the integer part of 1 € R). The time increments are chosen 
conform with the updating conventions of DSC schemes (as will be seen in a 
moment) and are consistent. In fact, in the first order of the time increment 
r and of the linear cell extension, the vector t V B Z in the centre point 
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of face l approximates the scalar products of the node vectors with the 
gradient VZ . More precisely, let for a fixed centre point on face t and 
e G M + the e-scaled cell have edge vectors t e~ : = e t e . Let also ,,V B Z M 
denote function (fTi)) for the e-scaled cell (with node vectors ^b^ = e^b). 
Then at the fixed point holds 

(15) < AI 6,grad(Z) > = „b ■ VZ = lim lim ^ t V s ~X , 

e—»0 r— »o e 

as immediately follows from the required C 1 -smoothness of the field Z. 

To recover the gradient VZ from (fl4|) in the same order of approxima- 
tion, observe that for every orthonormal basis ( v u) v =a,...,m-i of M m or C m , 
and for any basis (^6)^=0,. with coordinate matrix = < u u , ^b >, 
the scalar products of every vector a with ^6 equal 

(16) < M 6, a > = < „u > < „u, a > = /3£ a v 





( at the right-hand side, and henceforth, we observe Einstein's summation 
convention - yet without summing over indices that anywhere appear also 
as a left-hand subscript ) , hence 

(17) a v = tfa*, with ( 7 £) := ((^)*) _1 . 

This applied to the node vector basis ^b and (fT5|) yields the approximate 
gradient of Z at face t 

(18) ,VZ, = it ,V B ^- 

The scalar product of the gradient with face vector L f v = < J, „u > , 
v e {0, 1,2} is thus 

(19) t S = J ■ NZ = J v it t V B Z M = ^ N B Z^ . 

= : <,sf 

Continuity of the gradient at cell interfaces yields linear updating equations 
for Z p on the two adjacent faces. In fact, for any two neighbouring cells C, 
X with common face, labelled i in cell £ and n in x , continuity requires 

(20) lS = -lS. 

Substituting ([TO)) for t 5 and *5 and observing the time shifts in Q14p 
provides the updating relations for Z p at the cell interfaces. To make these 
explicit, we first introduce the quantities L zE ,n , ( t = 0, 5 ; /i = 0, 1, 2 ) 

' 2{-\y Z n ]t if n = [l/2] 

(2^+iZ'p - 2 „ZP) lt _ T/2 else 



(21) a n At) : = 



which in virtue of ([TJ) yields L z p — (p , Z) = Z o ,cr _1 ) = 

= t z^ 1 | Z o L o~ x , where (U : u i-> p denotes the nodal shift pertinent to 

face i . In particular 

(22) al m (*) = 2(-l) t ^| i , 

which together with (|2"Tj) is consistent for /i 7^ [t/2] with 
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(23) ^(t + r/2) = -\{2» +1 z? L + 2^){t). 

From (El [HI El [22]) follows that 

,S\ t+ r = ^{ L z^ l]t+T/2 - 2(-l)^/ 2 l L Z* ]t+T ) 

= ^(^|t + r/2-^ /2I ^l*+r)- 



(24) 



Continuity of Z , i.e. \z p = X K ZP , with HB) then implies 



(25) :< /21 (t+r) = 



*">2(* + t/2) + "a" >JJ(t + r/2) 



/2] ^ > « S ['/2] + (-!)<.+« * S [ K / 2 1 



For completeness we agree upon setting \z v ( i + r ) : = ^ z ™ ( £ + r/2 ) for 
/i ^ [t/2] (although this contains a slight inconsistency in that continuity 
might be infringed; this can easily be remedied by taking the arithmetic 
mean of the two adjacent values). - In fact, our agreement doesn't do harm, 
since any discontinuity disappears with mesh refinement. 

We dispose, hence, of a complete set of recurrence relations for z p ( given 
z n by the former reflection step) which at the same time determine the field 
components on face t and their gradients 

(26) NZ V = Us*. 

Essentially this constitutes the connection step of the algorithm. 
Nodal gradients are similarly (and even more simply) derived, using 

V B Z;(i+^) := ( W Z* - 2^){t); M = 0,1,2 

in the place of (fl4|) , and then again (fl8|) . With the node and cell-boundary 
values and gradients of T and u the nodal updating relations for the 
latter are immediately extracted from equations ([Til EH) m section [4]. For 
equation (|11|) this is essentially (up to the convective term) carried out in 
[Helj . section 5, and the procedure remains straightforward in the case at 
hand. Note that a well-timed LES coarsening routine, cf. section [3l should 
be periodically carried out before the nodal step of the iteration cycle in 
order to obviate instablities called forth by the energy cascade jPoj . 

The updating relations thus obtained are explicit and consistent with 
near-field interaction (only adjacent quantities enter). So, they can op- 
tionally be transformed into scattering relations for incident and reflected 
quantities ([6]) along the guidelines of section[2]- with established advantages 
for stability estimates [He2j . 



6. Pressure 



Unlike compressible flow - wherein a thermodynamic state equation relates 
pressure to density and temperature - Boussinesq-incompressible flow is 
conserved by pressure acting like a potential against violations of the local 
flow balance. Hence, the pressure is directly coupled to the flow divergence 
by Poisson's equation: 
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(27) 



= {qoo/t) diva.. 



The solution p for fluid velocities updated in the connection step provides a 
pressure gradient, which in the next reflection step - via equations (|10I12[) - 
repairs sporadic violations of the flow balance. Systematic test computations 
have shown that additional divergence clearing, as proposed by the author 
in earlier papers, e.g. [He3], is not necessary (nor even profitable) - and may 
be rather on the debit of computational performance and stability. 

In integral form, after using Gauss' Theorem, equation (|27[) becomes 



which can be solved by a Gauss-Seidel routine or by successive overrelax- 
ation, carried out after updating the face velocities in the connection step. 

In each iteration, firstly the discrete cell boundary integral Iq^ := 
J d( , u • dF is computed and explicitely - by solving a linear equation in 
Ig^ and the cellular pressures - the unique p n that exactly compensates 
Iq^ so that (|28|) is satisfied in every cell. In a second run over the mesh, 
continuity of the pressure gradient is re-established at the cell interfaces by 
updating the face pressures along the lines of the last section (viz. taking 
Z as p there). The loop is reiterated until (usually after a few iterations) 
equations (f2"5| hold with sufficient precision. 

7. Conclusion and Completions 

The DSC approach to computational fluid dynamics allows for simulation of 
turbulent flow in a plain and natural way. Cellular coarse-graining efficiently 
prevents turbulent pile up and preserves the large eddies. 

Since the energy cascade starts long before transition to turbulence 
(actually it is effective over the entire energy spectrum [Poj . [BIL1 pp. 72 ff.] ), 
a well-timed coarsening routine should escort every DSC fluid flow algo- 
rithm. Periodic coarse-graining significantly improves algorithm stability, 
even in the laminar regime - this without overly altering the flow dynamics, 
if the coarsening period is chosen as outlined in section [3] Handy criteria 
for a good choice - in relation to mesh refinement, e.g. - should be subject 
to further study. 

Details of implementation and special applications are clearly beyond 
the scope of the present study. However, we leave off with some examples 
for the purpose of illustration. 

The graphics displayed have been computed with author's test program 
DANSE. The latter combines a TLM Maxwell field solver with a DSC heat 
transfer and fluid flow algorithm written in the lines of this paper. 

In particular, the computations for coaxial RF power transmission line 
RL100-230 (motivated by an ion cyclotron resonance heating experiment 
in plasma physics) are in excellent keeping with reference computations 
(Fluent) and empirical data: For a line with inner conductor made of cop- 
per and air dielectric at atmospheric pressure the heating process has been 



(28) 
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simulated in horizontal position from standby to steady state operation with 
160 KW transmitted CW power at 100 MHz frequency, with outer conduc- 
tor cooled at 40 degrees Celsius temperature. Simultaneously, a Maxwell 
field TLM algorithm run in the same mesh provided the skin effect heat 
sources. 




Fig. 2. Jet over backward facing step 
{ highly turbulent regime ]. 




Fig. 3. Karman vortex street behind cylinder 
{ snapshot of oscillating flow ]. 
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Fig. 4. Air convection in horizontal coaxial line RL100-230 ; 
transverse section [ reference arrow: 0.1 ms~ ; transmitted power 160 kW CW, 
frequency 100 MHz; air dielectric; inner conductor copper, outer conductor cooled 
at 40 degrees Celsius ]. 
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Fig. 5. Heating process in coaxial line RL100-230 ; 
temperature maximum vs. time from power-on, 
convection compared to - fictitious - pure conduction 
[ transmitted power 160 kW CW, frequency 100 MHz, 
outer conductor cooled at 40 degrees Celsius ]. 
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Fig. 6. Cross sectional temperature profile in coaxial line RL100-230 computed at 
steady state [ transmitted power 160 kW CW, frequency 100 MHz; air dielectric; 
inner conductor copper, outer conductor cooled at 40 degrees Celsius ]. 



